#include "meshTypes.hpp"
#include "meshTypeFunctions.hpp"
#include "meshFunctions.hpp"
#include "structMeshFunctions.hpp"
#include "field.hpp"
#include "constantField.hpp"
#include "simulationTime.hpp"
#include "bc_base_all.hpp"
#include "timer.hpp"

#include "testHelpers.hpp"

#include <catch2/catch_test_macros.hpp>
#include <iostream>
#include <memory>

using namespace hms;

TEST_CASE( "StaticValue" ){

	auto test = []( auto&& dummyMesh, Index nComponents ){

		using MeshType = remove_qualifiers<decltype(dummyMesh)>;
		using BCType = StaticValue<MeshType>;

		auto mesh { uniMeshAs<MeshType>(3, 4, 1, 1, true) };

		const BoundaryPatch& patch { mesh.boundaryPatch("left") };

		Field f { mesh, nComponents };
		f.values() = 99;

		SECTION( "Uniform value"
			", nComponents=" + std::to_string(nComponents)
		){
			ArrayXs value (nComponents);
			value.setLinSpaced(nComponents, 1, nComponents);

			BCType& s {
				f.template setBoundaryCondition<BCType>(patch, value, true)
			};

			REQUIRE( s.isUniform() );

			for ( const auto& col :
				f.ghostValues(patch).colwise()
			)
				REQUIRE( (col == value).all() );

			REQUIRE( (s.get() == value).all() );
			REQUIRE( (s.getUniform() == value).all() );
		}
		SECTION( "List of values"
			", nComponents=" + std::to_string(nComponents)
		){
			ArrayXXs values ( nComponents, patch.size() );
			values.setRandom();

			BCType& s {
				f.setBoundaryCondition( BCType{f, patch, values, false} )
			};

			REQUIRE( !s.isUniform() );

			REQUIRE( (f.ghostValues(patch) == values).all() );

			REQUIRE( (s.get() == values).all() );
			REQUIRE( (s.getList() == values).all() );
		}
	};
	auto testSizes = [&]( auto&& dummyMesh){
		for (int nc=1; nc<=2; ++nc){
			test(std::move(dummyMesh), nc);
		}
	};
	TEST_ALL_MESHTYPES(testSizes)
}

TEST_CASE( "ZeroGradient" ){

	auto test = []( auto&& dummyMesh, Index nComponents ){

		using MeshType = remove_qualifiers<decltype(dummyMesh)>;
		using BCType = ZeroGradient<MeshType>;

		auto mesh { uniMeshAs<MeshType>(3, 4, 1, 1, true) };
		
		const BoundaryPatch& patch { mesh.boundaryPatch("left") };

		Field f { mesh, nComponents };

		BCType& z { f.template setBoundaryCondition<BCType>(patch) };

		auto check = [&](){
			REQUIRE( ( f.ghostValues(patch) == f.domainValues(patch) ).all() );
			/* for multiple ghost cell layers */
			// for ( auto row : patch.ghostCells().rowwise() )
			// 	REQUIRE(
			// 		(f.values()(hms::all, row) ==
			// 		f.domainValues(patch) ).all() );
		};

		f.values().setRandom();
		z.update();
		check();

		f.values().setRandom();
		z.update();
		check();
	};
	auto testSizes = [&]( auto&& dummyMesh ){
		for (int nc=1; nc<=2; ++nc){
			test(std::move(dummyMesh), nc);
		}
	};
	TEST_ALL_MESHTYPES(testSizes)
}

TEST_CASE( "DynamicValue" ){

	auto test = []( auto&& dummyMesh, Index nComponents ){

		using MeshType = remove_qualifiers<decltype(dummyMesh)>;
		using BCType = DynamicValue<MeshType>;

		auto mesh { uniMeshAs<MeshType>(3, 4, 1, 1, true) };
		
		const BoundaryPatch& patch { mesh.boundaryPatch("left") };

		Field f { mesh, nComponents };
		f.values() = 2;

		Field z { mesh, 1 };
		z.values() = 1;

		Time t { 10, 10 };

		BCType& d { f.template setBoundaryCondition<BCType>(patch, t) };

		SECTION( "Uniform time series"
			", nComponents=" + std::to_string(nComponents)
		){
			/* time dependent values */
			TimeSeries uniformVals (3);
			scalar timeIndex {0};
			
			for ( auto& pair : uniformVals ){
				pair.first = timeIndex;
				pair.second.resize(nComponents, 1);
				/* new values every two seconds, set to that time -5 */
				pair.second.setConstant( timeIndex - 5 );
				timeIndex += 2;
			}
			d.values( std::move(uniformVals), true );

			d.update();
			
			REQUIRE( (f.ghostValues(patch) == -5 ).all() );

			t.advanceBy( 3 );
			d.update();
			
			REQUIRE( (f.ghostValues(patch) == -3 ).all() );

			t.advanceBy( 1.5 );
			d.update();
			
			REQUIRE( (f.ghostValues(patch) == -1 ).all() );
		}
		SECTION( "Non-uniform time series"
			", nComponents=" + std::to_string(nComponents)
		){
			/* time dependent values */
			TimeSeries nonUniformVals (3);
			scalar timeIndex {0};
			
			for ( auto& pair : nonUniformVals ){
				pair.first = timeIndex;
				pair.second.resize( nComponents, patch.ghostCells().size() );
				/* new values every two seconds, set to that time -5 */
				pair.second.setRandom();
				timeIndex += 2;
			}
			d.values( std::move(nonUniformVals), false );
			d.update();

			REQUIRE( (f.ghostValues(patch) == d.values()[0].second).all() );

			t.advanceBy(2.5);
			d.update();

			REQUIRE( (f.ghostValues(patch) == d.values()[1].second).all() );

			t.advanceBy(2.5);
			d.update();

			REQUIRE( (f.ghostValues(patch) == d.values()[2].second).all() );
		}
		SECTION( "Uniform function"
			", nComponents=" + std::to_string(nComponents)
		){
			std::vector<FunctionParserWrapper> fparsers (nComponents);
			ArrayXs vals (nComponents);

			for ( auto& parser : fparsers ){
				parser.parse( "sin(t) + dt" );
			}
			d.functions( std::move(fparsers), true );

			auto req = [&](){
				vals = sin( t.current() ) + t.step();
				for ( auto col : f.values()
					(hms::all, patch.ghostCells().reshaped() ).colwise()
				)
					REQUIRE( col.isApprox(vals) );
			};

			d.update();
			req();
			t.advanceBy(2);
			d.update();
			req();
			t.advanceBy(1.4);
			d.update();
			req();
		}
		SECTION( "Non-uniform function"
			", nComponents=" + std::to_string(nComponents)
		){
			std::vector<FunctionParserWrapper> fparsers (nComponents);
			
			for ( auto& parser : fparsers ){
				parser.parse(
					"x*y + z + t"
					"- (x_dom+y_dom) + 2*z_dom"
					"+ val_dom*dt" );
			}
			d.functions( std::move(fparsers), false );
			d.zField(z);

			auto req = [&](){
				for ( Index g : patch.ghostCells() ){
					auto gc { mesh.cell(g) };
					auto dc { gc.domainCell() };
					ArrayXs vals (nComponents);
					vals =
						gc.centroid().prod() + z.values()(0, g) + t.current()
						- dc.centroid().sum() + 2*z.values()( 0, dc.index() )
						+ f.values()(0, dc.index() ) * t.step();

					REQUIRE( f.values().col(g).isApprox(vals) );
				}
			};

			d.update();
			req();
			t.advanceBy(2);
			d.update();
			req();
			t.advanceBy(1.5);
			d.update();
			req();
		}
	};
	auto testSizes = [&]( auto&& dummyMesh ){
		for (int nc=1; nc<=2; ++nc){
			test(std::move(dummyMesh), nc);
		}
	};
	TEST_ALL_MESHTYPES(testSizes)
}


TEST_CASE( "std::vector of boundary conditions" ){
	using MeshType = UniMesh;
	using BCType = BoundaryCondition<MeshType>;
	std::vector<std::unique_ptr<BCType>> b;
	Field<MeshType> f;
	BoundaryPatch p;
	Time t;
	b.push_back( std::make_unique<StaticValue<MeshType>>(f, p, true) );
	b.push_back( std::make_unique<ZeroGradient<MeshType>>(f, p) );
	b.push_back( std::make_unique<DynamicValue<MeshType>>(f, p, t) );
	REQUIRE( b[0]->type() == boundaryCondition::Type::staticValue );
	REQUIRE( b[1]->type() == boundaryCondition::Type::zeroGradient );
	REQUIRE( b[2]->type() == boundaryCondition::Type::dynamicValue );
}

TEST_CASE( "Adding boundary conditions to field" ){

	auto test = [&]( auto&& dummyMesh ){

		using MeshType = remove_qualifiers<decltype(dummyMesh)>;

		auto mesh { uniMeshAs<MeshType>(3, 4, 1, 1, true) };

		// Index ng { mesh.nGhostLayers() };
		Field field { mesh, 2 };
		field.values() = signalCode::uninitialised;
		field.domainValues() = 10;

		const BoundaryPatch
			& bottom { mesh.boundaryPatch("bottom") },
			& left   { mesh.boundaryPatch("left"  ) },
			& right  { mesh.boundaryPatch("right" ) },
			& top    { mesh.boundaryPatch("top"   ) };

		/****************/
		/* ZeroGradient */
		/****************/
		/* ZeroGradient does not require additional parameters */
		field
			.template setBoundaryCondition<ZeroGradient<MeshType>>(top)
			.update();

		/* must have the same value as domain field */
		REQUIRE( (field.ghostValues( top ) == 10).all() );

		/***************/
		/* StaticValue */
		/***************/
		Array2s uni { 1.5, 2.5 };
		field.template setBoundaryCondition<StaticValue<MeshType>>(
			bottom, uni, true
		);

		for ( const auto& col : field.ghostValues(bottom).colwise() )
			REQUIRE( (col == uni).all() );

		Array2Xs nonUni (2, 4);
		nonUni <<
			-1.5, math::e,  phys::g, 9.0,
			99.9, math::pi, 1234423, 42;

		field.template setBoundaryCondition<StaticValue<MeshType>>(
			right, nonUni, false
		);
		REQUIRE( (field.ghostValues(right) == nonUni).all() );

		/****************/
		/* DynamicValue */
		/****************/
		/* start at t = 2s, so that a time-dependent function can be tested */
		Time t { 2, 10, 10, scalar{10} };
		Field z { mesh, 1 };
		z.values() = 3;

		/* time-dependent values */
		TimeSeries
		uniT { { 0, Array2s{ 1234, 2345 } }, { 5, Array2s{ 6789, 7890 } } },
		uniTCopy {uniT},
		nonUniT { { 0, nonUni*2 }, { 5, nonUni*3 } },
		nonUniTCopy {nonUniT};

		field.template setBoundaryCondition<DynamicValue<MeshType>>(
			left, t, std::move(uniTCopy), true
		).update();

		for ( const auto& col : field.ghostValues(left).colwise() )
			REQUIRE( (col == uniT[0].second).all() );

		field.template setBoundaryCondition<DynamicValue<MeshType>>(
			right, t, std::move(nonUniTCopy), false
		).update();

		REQUIRE( (field.ghostValues(right) == nonUni*2).all() );

		/* functions */
		std::vector<FunctionParserWrapper>
			uniFunc { { "sin(t)" }, { "cos(t)"} },
			nonUniFunc { { "if(x*y<2, 0, 1)" }, { "z" } };

		/* uniform values set by function */
		field.template setBoundaryCondition<DynamicValue<MeshType>>(
			bottom, t, std::move(uniFunc), true, z
		).update();
		REQUIRE(
			( field.ghostValues(bottom).row(0) == sin(t.current()) ).all()
		);
		REQUIRE(
			( field.ghostValues(bottom).row(1) == cos(t.current()) ).all()
		);

		/* non-uniform values set by function */
		field.template setBoundaryCondition<DynamicValue<MeshType>>(
			top, t, std::move(nonUniFunc), false, z
		).update();
		for ( int i{0}; i<top.ghostCells().size(); ++i ){
			const auto& col { field.ghostValues(top).col(i) };
			const auto& cell { mesh.cell( top.ghostCells().reshaped()(i) ) };
			REQUIRE( col.x() == ( cell.centroid().prod() < 2 ? 0 : 1 ) );
			REQUIRE( col.y() == 3 );
		}
	};
	TEST_ALL_MESHTYPES(test)
}
